library(foreign)
library(gregmisc)
library(apsrtable)
library(Zelig)
library(grunfeld)
library(systemfit)
library(AER)
library(sem)
library(cem)
library(MatchIt)
library(Amelia)
library(DAAG)
library(bootstrap)

testmerge9 <- read.csv("baggott-charnysh.csv")

#####################################################
# Creating one T var from binary T0, T1, T2, T3, T4 #
#####################################################

testmerge9$T <- NA
	  
for(i in 1:nrow(testmerge9)){
    if (is.na(testmerge9$T1[i])==TRUE)    
 	  {testmerge9$T[i] <- NA             
 	  }
 	 else if (testmerge9$T1[i]==1)          
 	  {testmerge9$T[i] <- 1               
 	  }
 	  }
for(i in 1:nrow(testmerge9)){

    if (is.na(testmerge9$T2[i])==TRUE)     
     	  {testmerge9$T[i] <- NA             
 	  }
 	 else if (testmerge9$T2[i]==1)          
 	  {testmerge9$T[i] <- 2          }
 	  }

for(i in 1:nrow(testmerge9)){

    if (is.na(testmerge9$T3[i])==TRUE)     
 	  {testmerge9$T[i] <- NA              
 	  }
 	 else if (testmerge9$T3[i]==1)          
 	  {testmerge9$T[i] <- 3               
 	  }
 	  }

for(i in 1:nrow(testmerge9)){

    if (is.na(testmerge9$T4[i])==TRUE)     
 	  {testmerge9$T[i] <- NA              
 	  }
 	 else if (testmerge9$T4[i]==1)          
 	  {testmerge9$T[i] <- 4               
 	  }
 	  }  

for(i in 1:nrow(testmerge9)){

    if (is.na(testmerge9$T0[i])==TRUE)    
 	  {testmerge9$T[i] <- NA              
 	  }
 	 else if (testmerge9$T0[i]==1)          
 	  {testmerge9$T[i] <- 0               
 	  }
 	  }
testmerge9$T 

# Note that some nations are sequentially elected, e.g. Panama in 1975.  Then, T4=1 and T1=0, and T=4.We overwrite such that T=0 in these cases.  Just put T0 loop last in the list.

#########################################
# Checking for outliers driving results #
#########################################

# Democracy - no obvious outliers
plot(testmerge9$delta4demaut)


# GDP - there are some significant outliers
plot(testmerge9$delta4GDPpcWB) 

# Dropping the outliers:
outlierdrop <- testmerge9[testmerge9$delta4GDPpcWB < 100 & testmerge9$delta4GDPpcWB > -50,]
plot(outlierdrop$delta4GDPpcWB ~ outlierdrop$year)

# Aid-GDP pernicious effect still holds after dropping the outliers

lm(delta4GDPpcWB ~ totaid96, data=outlierdrop)
lm(delta4GDPpcWB ~ totaid96, data=testmerge9) 
lm(delta4GDPpcWB ~ scmem, data=testmerge9)
lm(delta4GDPpcWB ~ scmem, data=outlierdrop)

# Test for year outliers:
yeardrop <- NULL
yeardrop <- testmerge9[testmerge9$year!=1988 & testmerge9$year!=1989 & testmerge9$year!=1990 & testmerge9$year!=1991 & testmerge9$year!=1992 & testmerge9$year!=1993 & testmerge9$year!=1994 & testmerge9$year!=1995,]
plot(yeardrop$delta4GDPpcWB ~ yeardrop$year)


############################
# Summary statistics #
############################


a <- matrix(c(
# DVs
mean(na.omit(testmerge9$delta4GDPpcWB)),
sd(na.omit(testmerge9$delta4GDPpcWB)),
length(na.omit(testmerge9$delta4GDPpcWB)),

mean(na.omit(testmerge9$delta4demaut)),
sd(na.omit(testmerge9$delta4demaut)),
length(na.omit(testmerge9$delta4demaut)),

# Xs
mean(na.omit(testmerge9$AID)),
sd(na.omit(testmerge9$AID)),
length(na.omit(testmerge9$AID)),

mean(na.omit(testmerge9$us.totaid)),
sd(na.omit(testmerge9$us.totaid)),
length(na.omit(testmerge9$us.totaid)),

mean(na.omit(testmerge9$ciaaid.ussr)),
sd(na.omit(testmerge9$ciaaid.ussr)),
length(na.omit(testmerge9$ciaaid.ussr)),

# Controls
mean(na.omit(testmerge9$GDPpcWB)),
sd(na.omit(testmerge9$GDPpcWB)),
length(na.omit(testmerge9$GDPpcWB)),

mean(na.omit(testmerge9$demaut)),
sd(na.omit(testmerge9$demaut)),
length(na.omit(testmerge9$demaut)),

mean(na.omit(testmerge9$bigwar)),
sd(na.omit(testmerge9$bigwar)),
length(na.omit(testmerge9$bigwar)),

# Instruments
mean(na.omit(testmerge9$lpop)),
sd(na.omit(testmerge9$lpop)),
length(na.omit(testmerge9$lpop)),

mean(na.omit(testmerge9$britcol)),
sd(na.omit(testmerge9$britcol)),
length(na.omit(testmerge9$britcol))

), ncol=3, byrow=TRUE)
colnames(a) <- c("Mean","Standard Deviation", "Obs.")
library(xtable)
xtable(a)

################
# CEM Matching #
################

# 1 - Creating high/low aid bins

testmerge9$GDPpcWB <- as.numeric(testmerge9$GDPpcWB)
AidGDP <- cbind(orig$ccode, orig$year, orig$AidGDP)
colnames(AidGDP) <- c("ccode", "year", "AidGDP")
testmerge9 <- merge(testmerge9, AidGDP, by=c("ccode", "year"), all.x=TRUE)

testmerge9$hiaid <- NULL

for(i in 1:nrow(testmerge9)){

    if (is.na(testmerge9$AidGDP[i])==TRUE)     # if GDP is NA
 	  {testmerge9$hiaid[i] <- NA}              # then hi is NA
 	 else if (testmerge9$AidGDP[i]>=7)          # if GDP>3000
 	  {testmerge9$hiaid[i] <- 1}                # then hi is 0
 	 else if (testmerge9$AidGDP[i]<7)          # if GDP>3000
 	  {testmerge9$hiaid[i] <- 0}  
 	  }

testmerge9$AID <- as.numeric(testmerge9$AID)
testmerge9$hitotaid <- NULL
for(i in 1:nrow(testmerge9)){
	if(is.na(testmerge9$AID[i])==TRUE)
	{testmerge9$hitotaid[i] <- NA}
	else if (testmerge9$AID[i] >=2500 )
	{testmerge9$hitotaid[i] <- 1}
	else if(testmerge9$AID[i] < 2500 )
	{testmerge9$hitotaid[i] <- 0}
	}
summary(testmerge9$AID)

testmerge9$us.totaid <- as.numeric(testmerge9$us.totaid)
testmerge9$hiusaid <- NULL
for(i in 1:nrow(testmerge9)){
	if(is.na(testmerge9$us.totaid[i])==TRUE)
	{testmerge9$hiusaid[i] <- NA}
	else if (testmerge9$us.totaid[i]>=500 )
	{testmerge9$hiusaid[i] <- 1}
	else if(testmerge9$us.totaid[i] < 500 )
	{testmerge9$hiusaid[i] <- 0}
	}

test <- as.data.frame(cbind(testmerge9$AidGDP, testmerge9$hiaid))
test[1:500,]
# for loop worked.

# 2 - get coarsened exact matching SATT
# Approach: Match based on starting values of GDP per capita and democracy; X is high aid flows and Y is change in GDPpc or democracy.

# 2.1 - Restriction to SC country-terms, X=us.totaid, Y= change GDP

temps <- as.data.frame(testmerge9[testmerge9$scmem==1,])

cemset <- cbind(temps$GDPpcWB, as.numeric(temps$hiusaid), temps$demaut, temps$delta4GDPpcWB, temps$ccode, temps$year, temps$lpop, temps$britcol)
colnames(cemset) <- c("GDPpcWB", "hiaid", "demaut", "delta4GDPpcWB", "ccode", "year", "lpop", "britcol")
head(cemset)
na.cemset <- as.data.frame(na.omit(cemset))
mat <- cem(treatment = "hiaid", data=na.cemset, drop=c("delta4GDPpcWB", "ccode", "year", "lpop", "britcol"))
mat
cem.match.att <- att(obj=mat, formula=delta4GDPpcWB ~ hiaid + GDPpcWB + demaut, data = na.cemset, model="linear")
cem.match.att
xtable(summary(cem.match.att))

pre.imbalance <- imbalance(group=na.cemset$hiaid, data=na.cemset, drop=c("hiaid", "ccode", "year", "lpop", "britcol", "delta4GDPpcWB"))
pre.imbalance

# 2.2 - Restriction to SC country-terms, X=AID, Y= change GDP

cemset <- cbind(temps$GDPpcWB, temps$hitotaid, temps$demaut, temps$delta4GDPpcWB, temps$ccode, temps$year, temps$lpop, temps$britcol)
colnames(cemset) <- c("GDPpcWB", "hiaid", "demaut", "delta4GDPpcWB", "ccode", "year", "lpop", "britcol")
head(cemset)
na.cemset <- as.data.frame(na.omit(cemset))
mat <- cem(treatment = "hiaid", data=na.cemset, drop=c("delta4GDPpcWB", "ccode", "year", "lpop", "britcol"))
mat
cem.match.att <- att(obj=mat, formula=delta4GDPpcWB ~ hiaid + GDPpcWB + demaut, data = na.cemset, model="linear")
cem.match.att
xtable(summary(cem.match.att))

pre.imbalance <- imbalance(group=na.cemset$hiaid, data=na.cemset, drop=c("hiaid", "ccode", "year", "lpop", "britcol", "delta4GDPpcWB"))
pre.imbalance

# 2.3 - Restriction to SC country-terms, X=AID, Y= change dem

cemset <- cbind(temps$GDPpcWB, temps$hitotaid, temps$demaut, temps$delta4demaut, temps$ccode, temps$year, temps$lpop, temps$britcol)
colnames(cemset) <- c("GDPpcWB", "hiaid", "demaut", "delta4demaut", "ccode", "year", "lpop", "britcol")
head(cemset)
na.cemset <- as.data.frame(na.omit(cemset))
mat <- cem(treatment = "hiaid", data=na.cemset, drop=c("delta4demaut", "ccode", "year", "lpop", "britcol"))
mat
cem.match.att <- att(obj=mat, formula=delta4demaut ~ hiaid + GDPpcWB + demaut, data = na.cemset, model="linear")
cem.match.att
xtable(summary(cem.match.att))
# Total aid has positive effect sig at 10%.

pre.imbalance <- imbalance(group=na.cemset$hiaid, data=na.cemset, drop=c("hiaid", "ccode", "year", "lpop", "britcol", "delta4demaut"))
pre.imbalance

######## 2.4 - Restriction to SC country-terms, X=us.totaid, Y= change dem

cemset <- cbind(temps$GDPpcWB, temps$hiusaid, temps$demaut, temps$delta4demaut, temps$ccode, temps$year, temps$lpop, temps$britcol)
colnames(cemset) <- c("GDPpcWB", "hiaid", "demaut", "delta4demaut", "ccode", "year", "lpop", "britcol")
head(cemset)
na.cemset <- as.data.frame(na.omit(cemset))
mat <- cem(treatment = "hiaid", data=na.cemset, drop=c("delta4demaut", "ccode", "year", "lpop", "britcol"))
mat
cem.match.att <- att(obj=mat, formula=delta4demaut ~ hiaid + GDPpcWB + demaut, data = na.cemset, model="linear")
cem.match.att
xtable(summary(cem.match.att))
# US aid has no significant effect on democracy.
pre.imbalance <- imbalance(group=na.cemset$hiaid, data=na.cemset, drop=c("hiaid", "ccode", "year", "lpop", "britcol", "delta4demaut"))
pre.imbalance

#################################
# Plotting repeat members graph #
#################################


temps$startyear <- NA
for(i in 1:nrow(temps)){
	if(is.na(temps$T[i])==TRUE)
	{temps$startyear[i] <- NA}
	else if(temps$T[i]==1)
	{temps$startyear[i] <- as.vector(temps$year[i])}
	}

temps$endyear <- NA
for(i in 1:nrow(temps)){
	if(is.na(temps$T[i])==TRUE)
	{temps$endyear[i] <- NA}
	else if(temps$T[i]==2)
	{temps$endyear[i] <- as.vector(temps$year[i])}
	}

temps$endyearlag <- c(temps$endyear[2:length(temps$endyear)], NA)
matrix <- as.data.frame(cbind(as.vector(temps$ccode), temps$startyear, temps$endyearlag))

matrix.nona <- na.omit(matrix)
matrix.nona[,1] <- as.numeric(as.vector(matrix.nona[,1]))
matrix.nona[,2] <- as.numeric(as.vector(matrix.nona[,2]))
matrix.nona[,3] <- as.numeric(as.vector(matrix.nona[,3]))

count.country <- matrix(data=NA, ncol=2, nrow=length(unique(matrix.nona[,1])))
count.country[,1] <- unique(matrix.nona[,1])

for(i in 1:nrow(count.country)){
	count.country[i,2] <- length(matrix.nona[count.country[i,1]==matrix.nona[,1],1])
	}

unique.country.index <- count.country[count.country[,2]>1,]

unique.country.mat <- matrix.nona[matrix.nona[,1] %in% unique.country.index[,1],]

index.plot.country <- matrix(data=NA, ncol=2, nrow=length(unique(unique.country.mat[,1])))
index.plot.country[,2] <- unique(unique.country.mat[,1])
index.plot.country[,1] <- 1:nrow(index.plot.country)

full.plot.index <- rep(NA, nrow(unique.country.mat)) 
for(i in 1:nrow(unique.country.mat)){
	full.plot.index[i] <- index.plot.country[index.plot.country[,2]==unique.country.mat[i,1],1]
	}

pdf(file="repeats.pdf")
plot(x=2005, y=nrow(unique.country.index),  xlim=c(1946,2005), cex=0, ylim=c(0,nrow(unique.country.index)), yaxt="n",  main="Repeat Temporary Members", xlab="Year", ylab="Country")
segments(x0=matrix.nona[,2], y0=full.plot.index, x1=matrix.nona[,3], y1=full.plot.index, col="red")
dev.off()

# 66 countries served more than one term on the council.
# Only 10 cases in which a country served sequentially.

##########################
# Imputation with Amelia #
##########################

amelia9 <- cbind(
  as.numeric(testmerge9$year),
  as.numeric(testmerge9$ccode),
  as.numeric(testmerge9$unmem),
  as.numeric(testmerge9$scmem),
  as.numeric(testmerge9$democ),
  as.numeric(testmerge9$autoc),
  as.numeric(testmerge9$polity2),
  as.numeric(testmerge9$lpop),
  as.numeric(testmerge9$GDP2000dollarsWB),
  as.numeric(testmerge9$GDPpcWB),
  as.numeric(testmerge9$TI_corruption),
  as.numeric(testmerge9$totaid96),
  as.numeric(testmerge9$bigwar),
  as.numeric(testmerge9$us.totaid),
  as.numeric(testmerge9$AID),
  as.numeric(testmerge9$ciaaid.ussr),
  as.numeric(testmerge9$demaut),
  as.numeric(testmerge9$britcol)
)
colnames(amelia9) <- c("year", "ccode", "unmem", "scmem", "democ", "autoc", "polity2", "lpop", "GDP2000dollarsWB", "GDPpcWB", "TI_corruption", "totaid96", "bigwar", "us.totaid", "AID", "ciaaid.ussr", "demaut", "britcol")

a.out <- amelia(x=amelia9, cs="ccode", ts="year", m=10)
a.out
save(a.out, file = "imputations.RData")
write.amelia(obj=a.out, file.stem="imputations", format="dta")

# Each imputation takes one minute.
# Diagnostic tests:
plot(a.out, which.vars=9:18)
overimpute(a.out, var=16)  # imputing USSR aid looks risky...
overimpute(a.out, var=15) # don't know how to interpret this
missmap(a.out)
disperse(a.out)


write.amelia(obj=a.out, file.stem="outdata")

a1 <- read.csv("outdata1.csv")
a2 <- read.csv("outdata2.csv")
a3 <- read.csv("outdata3.csv")
a4 <- read.csv("outdata4.csv")
a5 <- read.csv("outdata5.csv")
a6 <- read.csv("outdata6.csv")
a7 <- read.csv("outdata7.csv")
a8 <- read.csv("outdata8.csv")
a9 <- read.csv("outdata9.csv")
a10 <- read.csv("outdata10.csv")

avg <- (a1+a2+a3+a4+a5+a6+a7+a8+a9+a10)/10
names(avg)
write.csv(avg, "imputationdata.csv")

# Exported average of 10 datasets.  Then re-ran IV regressions in stata.

##########################################
# Re-running regressions on imputed data #
##########################################

# Robust in stata.

##########################################
# Instrumental Variables Analysis        #
##########################################

m <- lm(lpop ~ us.totaid, data=testmerge9)
m <- lm(lpop ~ AID, data=testmerge9)
m <- lm(lpop ~ ciaaid.ussr, data=testmerge9)
m <- lm(lpop~ ciaaid.prc, data=testmerge9)

m <- lm(britcol ~ us.totaid, data=testmerge9)
m <- lm(britcol ~ AID, data=testmerge9)
m <- lm(britcol ~ ciaaid.ussr, data=testmerge9)
m <- lm(britcol ~ ciaaid.prc, data=testmerge9)

m <- (delta4GDPpcWB ~ lpop, data=testmerge9)
m <- (delta4demaut ~ lpop, data=testmerge9)
m <- (delta4GDPpcWB ~ britcol, data=testmerge9)
m <- (delta4demaut ~ britcol, data=testmerge9)

summary(m)

fml <- list("mu" = delta4GDPpcWB ~ AID + GDPpcWB + demaut + bigwar,  "inst" = lpop ~ AID)
m <- zelig(formula = fml, model = "twosls", data=testmerge9)
summary(m)

m <- ivreg(delta4GDPpcWB ~ AID + GDPpcWB + demaut + bigwar | lpop, data=ivset)

m <- tsls(delta4GDPpcWB ~ AID + GDPpcWB + demaut + bigwar, instruments = ~ lpop, data=testmerge9)

ivset <- as.data.frame(cbind(testmerge9$delta4GDPpcWB, testmerge9$AID, testmerge9$GDPpcWB, testmerge9$demaut, testmerge9$bigwar, testmerge9$lpop))
colnames(ivset) <- c("delta4GDPpcWB", "AID", "GDPpcWB", "demaut", "bigwar", "lpop")


###########################
# K-Fold Cross Validation #
###########################

cv.lm(df=testmerge9, form.lm=formula(delta4demaut ~ hitotaid + GDPpcWB + demaut), m=3, dots=FALSE, seed=02138, plotit=TRUE, printit=TRUE)
# in DAAG package
# can also use crossval() in bootstrap package

